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Abstract 

The one-component plasma (OCP) represents the simplest statistical mechanical 
model of a Coulomb system. For this reason, it has been extensively studied over 
the last forty years. The advent of the integral equations has resulted in a dramatic 
improvement in our ability to carry out numerical calculations, but came at the 
expense of a physical insight gained in a simpler analytic theory. In this paper we 
present an extension of the Debye-Hiickel (DH) theory to the OCP. The theory 
allows for analytic calculations of all the thermodynamic functions, as well as the 
structure factor. The theory explicitly satisfies the Stillinger-Lovett and, for small 
couplings, the compressibility sum rules, implying its internal self consistency. 

Key words: One-component plasma, electrolytes, structure factors, Debye-Hiickel 
theory. 



1 Introduction 

The classical one-component plasma (OCP) is an idealized system of N identical point-particles 
of charge q, in a uniform neutralizing background of dielectric constant D [1,2]. For concreteness 
we shall suppose that the particles are positively charged, while the background is negative. 
Each ion, inside the volume V, is assumed to interact with the others exclusively through the 
Coulomb potential. The OCP has been extensively studied because it serves as the simplest 
possible model for a variety of important physical systems, ranging from electrolytes and charge- 
stabilized colloids [3] to the dense stellar matter [4]. With the advent of powerful computers and 
new developments in the liquid-state theory our ability to perform thermodynamic calculations 
on this model has seen a dramatic improvement characterized, in particular, by the quantitative 
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agreement between the integral-equations based theories of the OCP and the Monte Carlo 
(MC) simulations [5] . Unfortunately, the intrinsic complexity of the integral equations, instead 
of clarifying the underlying physics, tends to obscure it. This should be contrasted with the 
simplicity and the transparency of the Debye-Hiickel (DH) theory [6], which allows for a very 
clear physical picture of an ionic solution. 

Let us recall the early history of electrolyte solutions. The first modern ideas about electrolytes 
can be traced to the pioneering work of the Swedish physicist Svante Arrhenius [7] at the end 
of the last century. In particular, it was Arrhenius who realized that when salts and acids are 
dissolved in a polar solvent, the molecules become dissociated, producing cations and anions. 
Arguing from what can now be called mean-field point of view, Arrhenius concluded that, since 
the electrolyte solution is charge neutral and the ions are uniformly distributed, the average force 
acting on each particle is null. All the nontrivial characteristics of the electrolytes Arrhenius 
attributed to the incomplete dissociation of the molecules. In this simple picture an electrolyte 
is treated as an ideal gas composed of three species, cations, anions, and neutral molecules, 
whose densities are controlled by the law of mass action. All the electrostatic interactions are 
neglected except in as far as treating the cations and the anions as distinct entities. This simple 
theory has proven to work quite well for what are now known as weak electrolytes, such as 
Br0nsted acids and bases. In the case of strong electrolytes, such as NaCl or HC1, the theory 
proved seriously flawed. It took almost forty years before a satisfactory solution could be found. 
It appeared in the form of the now famous DH theory of strong electrolytes [6] . The great insight 
of Debye and Hiickel was to realize that although the ions are on average uniformly distributed 
inside the solution, due to the long-range Coulomb force there exist strong correlations in the 
positions of the positively and the negatively charged particles; evidently in the vicinity of a 
positive ion there will be an excess of negative particles and vice-versa. 

To make this idea concrete and to see how it applies to the OCP, let us fix one mobile ion at 
the origin and ask what is the induced electrostatic potential in its surrounding. Clearly this 
must satisfy the Poisson equation, 

Att 

V 2 0(r) = --0(r). (1) 

To find the closure to this equation, we shall follow DH and suppose that the rest of the mobile 
ions arrange themselves in accordance with the Boltzmann distribution, 

e {r) = qp + exp [-Pq4>(r)\ - qp-, (2) 

where p+ = N/V is the average density of the mobile ions, p_ is the density of the uniform 
neutralizing background, and f3 = Xjk^T. The next step of the DH theory is to linearize the 
exponential factor, leading to 

Q(r) = -^W), (3) 
where k d = v /47tAbP+ and Ab = fiq 2 / D are the inverse Debye screening and the Bjerrum 
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lengths, respectively. Clearly, the linearization can only be justified in the high-temperature 
(weak-coupling) limit. The resulting Helmholtz equation, V 2 0(r) = K^<j)(r), can be easily 
integrated to produce a potential of the Yukawa form. The fundamental lesson of DH is that 
the ions arrange themselves in such a way as to screen the long-range Coulomb interaction. It 
is this renormalization of the interaction potential that is responsible for the existence of the 
thermodynamic limit for Coulomb systems. However, not everything is rosy with this simple 
theory. It is sufficient to look at the charge- density distribution, 



to notice that something went seriously wrong. 

Clearly, the physical restriction that g(r) > —qp- is strongly violated in the region near the 
fixed ion. The problem can be traced back to the linearization of the Boltzmann factor, which 
is unjustified at short distances, since there the electrostatic potential is not small even for 
high temperatures. Fortunately, not all is lost. A simple solution to overcome this difficulty 
was suggested by Nordholm [8], who proposed an augmentation of the DH theory to include 
an effective spherical cavity around the fixed ion, inside which no other ions can penetrate. 
The presence of such a cavity is quite reasonable, since the electrostatic repulsion between the 
like-charged ions should prevent them from coming into close contact. Furthermore, we can 
estimate the size of the hole, h, by comparing the repulsive Coulomb energy with the kinetic 
thermal energy, q 2 /Dh ~ ksT, or h ~ q 2 jDhsT. Evidently, the higher the temperature, the 
smaller the size of the exclusion region. This, of course, is intuitive, since at higher temperatures 
the ions will have more kinetic energy to overcome the mutual repulsion. A more consistent 
way of defining the hole size h is to require an overall charge neutrality [8], or, equivalently, 
the continuity of the potential and of the electric field across the hole boundary. Performing a 
simple calculation, we find 



/i = d[w(r)-i]/V3r, (5) 
cu(r) = [i + (3r) 3/2 ] 1/3 , (6) 

where we have defined the usual coupling constant for the OCP, V = (3q 2 /Dd, and d is the 

characteristic length scale, d = (Anp^/?^ . In the low-coupling limit this reduces to the 
energetically determined expression for the size of the cavity, h ~ Ab- Nordholm was able to 
demonstrate that this Debye-Hiickel plus hole (DHH) theory produces an equation of state for 
the OCP which is in good agreement with the MC simulations. The question, however, still 
remains to what extent the hole is a physical object or just a convenient mathematical trick 
to correct for the linearization of the Poisson-Boltzmann (PB) equation. Clearly, if the cavity 
postulated by the DHH theory is real, the best way to study it is by considering the structure 
factor. In particular, if everything is all right with the DHH theory, the structure factor obtained 
on its basis should be in good agreement with the MC simulations. Unfortunately, it is well 
known that the traditional ways of obtaining the correlation functions out of the DH theory 
lead to expressions which are seriously flawed [9]. In the case of the restricted primitive model 
(RPM), equisized spheres carrying charges ±q, the correlation functions violate the well-known 
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Stillinger-Lovett sum rules [10] and do not reproduce the charge-density oscillations known to 
exist at high densities [11]. 



Recently Lee and Fisher [9] have proposed an extension of the DH theory of strong elec- 
trolytes to nonuniform densities. The generalized DH theory (GDH) allows the calculation 
of the density-density and of the charge-charge correlation functions in a most natural way, 
through a functional differentiation of the free-energy functional. Furthermore, since the theory 
is constructed at the level of free energy, it is internally self consistent, as can be judged by 
the various sum rules that it satisfies. This, of course, is a great advantage over the traditional 
integral-equations based theories, which are constructed at the level of the correlation functions 
and depend on the route taken to the thermodynamics [12], i.e. virial, compressibility, etc. The 
comparison of the GDH theory with experiments or simulations, however, is made difficult by 
the same flaw (or virtue, depending on how one looks at it) that the original DH theory suffers, 
its linearity. The linearity of the DH theory for electrolytes is both a blessing and a curse. It 
ensures the internal self consistency of the theory, but is also responsible for undercounting the 
configurations in which the oppositely charged ions come into a close contact, forming dipolar 
pairs. It was shown recently how this difficulty can be overcome in the context of the Debye- 
Hiickel-Bjerrum (DHBj) theory [13], by allowing for the existence of a chemical equilibrium 
between the monopoles and dipoles. Unfortunately this stratagem is difficult to implement in 
the case of the GDH theory. The goal of this paper, then, is two-fold: test the physical nature 
of the cavity in the DHH theory and, by apply the GDH theory to the OCP — which is free 
from the cluster formations that plague RPM — test the extent of its validity. 



2 The generalized Debye-Hiickel theory for the one-component plasma 

The DHH theory for the OCP is extended to allow for a nonuniform, slowly varying ionic 
density, 



To preserve the electroneutrality on long-length scales, the overall equilibrium densities must 
be equal, p + = p_. 

The Helmholtz free-energy, F, is a functional of the ionic density p+{r). The direct correlation 
function, C ++ (ri — r 2 ), is defined in terms of the second derivative of the excess free energy, 



p+(r) = p + (1 + A cos k ■ r) . 



(7) 



The negative background, as in the original OCP theory, is maintained uniform, 



p_(r)=p_, Vr. 



(8) 



C ++ (n-r 2 ) = -/3 



{F [p + (r)} - F ideal [ P+ (r)]} 
8p+(ri)5p+(r 2 ) 



P+(r)=P+ 
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*(ri-r 2 ) p S^F[ P+ (r)] 



(9) 

P+(r)=P+ 



P+ Sp+(r l )5p + (r 2 ) 
Here, -Fldeai [p+(^)] is the usual ideal-gas Helmholtz free-energy functional, 

/SFideai [ P+ (r)\ = J d 3 r' P+ (r') {in [ P+ (r')A 3 ] - l} , (10) 

where A is the thermal de Broglie wavelength. 

The direct correlation function is connected with the total correlation function, H(r), through 
the Ornstein-Zernike relation, 

H(r) = C ++ (r) +p + j dV C ++ (r - r') H(r'), (11) 

which in the Fourier space can be written as 

H(k) = C+ ^ k) , (12) 



where C ++ (k) and H{k) are the Fourier transforms of the direct and the total correlation 
functions, respectively, 



C ++ (k) = J d 3 r C ++ (r) exp (ifc • r) , (13) 
H(k) = J d 3 r H(r) exp (ifc • r) . (14) 

The structure factor is defined as 

S(k) = 1 + p + H(k) = _ 1 (15) 

1 - P+C++{k) 

Evidently, the knowledge of the direct correlation function C ++ (r) is equivalent to the knowl- 
edge of the structure factor S{k). 

To proceed, we impose an infinitesimal variation on the mobile- ion density, Eq. (7), and expand 
the reduced Helmholtz free-energy functional density, / = (3F/V, in powers of A. To second 
order, the variation 5f can be written as (see details in appendix B) 

5f [ P+ (r)\ = f [ P+ (r)\ - f [p + ] = (3jlp + A5 k o + V x (fc)p + A 2 (1 + 5 k0 ) , (16) 

where f3p = df /dp + is the equilibrium chemical potential, S^o = -^(27r) 3 <5 3 (fc) is the Kronecker 
delta and <5 3 (fc) is the three-dimensional Dirac delta function. The free-energy density of the 
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homogeneous reference system, / [p+], is obtained by setting A = in the expression for 
/ [p+(r)]. Clearly, if we are able to construct the free-energy functional for a nonuniform system, 
the structure factor, S(k), can be read directly from the second-order term. 

We proceed in a way exactly analogous to the usual DH theory. Let us fix one positive ion at 
r' and ask what is the electrostatic potential, <f>(r, r'), at a point r in its surrounding. We shall 
assume that, just like in the uniform case, this potential satisfies the PB equation, 

V 2 Mr, ^ = ~ {* s (r - r') + p + (r) exp [-(3q^(r, r')} - p_(r)} . (17) 

A crucial point to remark [9] is that the potential which appears in the Boltzmann factor of (17), 
4>(r,r'), represents a local-induced potential, 

4>(r,r') = <j)(r,r') -$(r), (18) 



obtained by extracting from the total potential 4>(r, r') an imposed electrostatic potential, $(r), 
produced by the neutralizing background and the imposed charge-density variation (7), 

V 2 $(r) = -^[p + (r)-p_] = -^±Acosfc-r, Vr. (19) 

With the separation of the total potential (f>(r, r') into two parts, the Helmholtz free-energy 
functional can be written as 

F[p + (r)] = F idcal [p+(r)] + 

-^Imposed \p + (r)} + -^induced \p + (r)}, (20) 
where the excess free energies are obtained through the Debye charging process [6,9], 

i 

imposed [p + (r)\ = q J dV [p + (r') - p_] J d\ $(r', Xq), (21) 

o 

i 

induced \p+(r)] = q J dVp + (r') J d\^(r',Xq). (22) 



In Eq. (22), tp(r') is the mean induced electrostatic potential felt by the positive ion fixed at 



ip(r') = lim 

r— >r' 



4>(r, r') - 



D\r — r' 



(23) 



Linearization of the Boltzmann factor of (17) results in the GDH equation for the induced 
potential, 
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V 2 0(r, r') = \5\r - r') - p + {r) 



D 

4nq 

~D~ 



5 3 (r-r') -p+(l + Acosfc-r) , for \r - r'\ < h, (24) 



vli(r y )= ^^n kry) 

= Kn(l + A cos k ■ r)(f)(r,r'), for \r - r'\ > h. (25) 

As discussed in the introduction, to prevent the unphysical artifacts of the linearization of 
the PB equation, we have explicitly introduced a cavity h around the fixed ion at r', given 
by Eq. (5), into which no other mobile ions can penetrate. 

In the following subsections we shall obtain the contributions to the variation of the free-energy 
density. 



2.1 The ideal-gas contribution 

The ideal-gas contribution is given by (10) with the imposed mobile-ion charge distribution (7). 
Expanding (10) up to order A 2 and using the integrals (A. 3) and (A. 4), we obtain the ideal-gas 
contribution to the variation of the reduced free-energy density, 

5/idcai [p+{r)\ = In (p + A 3 )p + A4o + \p+A 2 (1 + 5 k0 ) . (26) 

2.2 The imposed electrostatic potential contribution 



The imposed electrostatic potential satisfies the Poisson equation (19), whose formal solution 
can be written as 

*W = >/dV^f (27) 

The contribution to the Helmholtz free-energy functional is obtained through the Debye charg- 
ing process [6], 



imposed [p + {r)\ = q J dV [p + (r') - p_] J d\ $(r', Xq) 



= A* / d 3 r d 3 r > COsk - r _ CO * k - r ' j dA A. (28) 



In this case the charging merely produces a trivial factor of 1/2 and using the integral (A. 7), we 
obtain the contribution of the imposed electrostatic potential to the variation of the reduced 
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free-energy density, 



^imposed [ P+ (r)} = \ ^^p + A 2 (l + 5 k0 ) 



(29) 



2.3 The induced electrostatic potential contribution 



The induced electrostatic potential satisfies the GDH equation, given by (24) and (25). It is 
convenient to rewrite them in a spherical coordinate system centered on the positive ion fixed 
at r'. Introducing the difference vector 



R = r - r', 



(30) 



the GDH equation for the induced potential reads 



, i-^{5"(R)-p + [l + Acosk.(R + r')}}, 

[I + A cos k-(R + r')] <j>{R + r', r'), 



for \R\ < h, 
for \R\> h. 



(31) 



Using the Green's function G (R, R') associated with (31) derived in appendix C, it can be 
transformed into the integral equation 



0(r, r') = 4>{R + r', r') = 1 J d 3 R> g (R f ) G (R, R') , 



(32) 



where the effective charge density, g(R), is given by 



q(R) 



qS 3 (R) - qp + [1 + A cos k ■ (R + r')] , for R < h, 



D 

Aix 



k14> (R + r', r') A cos k ■ (R + r'), for R>h. 



This equation can be solved perturbatively in powers of A (see appendix D). 
The mean induced electrostatic potential felt by the positive ion fixed at r*', 



ib(r') = lim 



(R + r'y) - 



q 



DR 



(33) 



(34) 



can be written, to order A 2 , as (see derivation in appendix D) 



f3qip(r') = —-x(x + 2) — A cos k ■ r 



1 sin ax cos ax x 

2 a 3 (l + x) a 2 (l + x) + l + x 



a 



H ( x i a ) 
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+ ~T^ A2 E (2* + 1) cos 2 (k-r' + e^) 3 ^+ l( ^ ) 2?(x, a) 
+ yf^A 2 £ (-if +1 (2€ + 1) cos 2 (k.r' + 



x [lgt+i(-x) 



T£(x, a) + a)l e (x, a) — 2%(x, a) \ . 



2 ge+i(x) 

where x = nj)h, a = k/nn, je(Q is the spherical Bessel function of the first kind, 



(35) 



7T 



(36) 



ge(£) is the £-th grade polynomial associated with the modified spherical Bessel function of the 
third kind, ke(£), 



m=0 



m=0 



2 m m! V 2m 



(37) 



where Y{m) = (m—l)\ is the Euler gamma function, and {X^}, v — ±, 0, are the one-dimensional 
quadratures 



o 

oo 

2^(x, a) = y ds s~ e g e (s)jt(as)le(s, a) exp [2 (a - s)] , 

X 

oo 

J/(x, a) = J ds s~ e ge{s)je(as) exp [2(z - s)] . 



(38) 
(39) 
(40) 



The contribution to the Helmholtz free-energy functional is obtained through the Debye charg- 
ing process [6,9], 



induced \p+(r)] = Q J dVp+(r') J dX^(r',Xq) 



1 

= QP+ J d 3 ^' (1 + A cos fc • r') j d\tp(r', Xq), 



(41) 



which yields the reduced free-energy density, 



/induced [P+( r )\ = P+fo + 



h + \fM) 



p + A5 k0 + l - [A (a) + f 2 (a)} p + A 2 (1 + 5 k0 ) , (42) 



from which the variation follows, 



5/induced \p+(r)] 



/induced [P+( r , 



— /induced [P+\ 

p+A5 k0 + \ [h{a) + f 2 (a)} p + A 2 (1 + 5 k0 ) , 



where 



u\ (u x + 1) 



^ 2 i -i 
L){ + L) X + 1 



1 

4 



2 2?r , ( W + W + 1 

1 — or H = + In 



3^ 



2 , (2uo+ 1 



1 2 r A 2 2 } A 

/i(a) = + — / dA — sin(a:r A / A) + — / dA — cos(ax A / A) 

cr a' 3 J ux a z J lu\ 



d\ u\ — 1 
A w A 



/ 2 (a) = E(^+ 1 ) 



1 + (-I)^m 



1 + 5 



fcO 



1 



- / dA — r — ll(x x ,a X) 



+ 



e+ i f dA lu\ - 1 



^+i(-^a) 
, ^+i(^a) 



+ 2I+(x x , a/X)li{xx, a/A) - 22£(x A , a/A) 



o-a = ^a - 1, 



1 + A 3 (3r) 



3/2 



1/3 



i + (3r) 3/2 



1/3 



We note that the reduced induced free-energy density for the reference system, 



/induced [p+] — P+fo 



2 2tt /cu 2 + cu + r 



2 (2uu + 1 N 



is the same as the one previously obtained by Nordholm [8]. 
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3 Analytical results and the sum rules 



Collecting all the contributions to the variation of the reduced free-energy density, 

S f [P+( r )} = 5 /ideal [P+(r)} + 5 /imposed [p+( r )] + ^/induced \p+(r)] , 



(51) 



and comparing with the expansion given by (16), we obtain the equilibrium chemical potential, 
Pfl = In (p + A 3 ) + f + \ lim h{a) = In (p+A 3 ) + f + 1 (l - uo 2 ) , (52) 



which corresponds to the usual OCP chemical potential [8], and the structure factor, 



2 r A 2 2 r A 

S ,-1 (fe) = 1 + — / dA — sin(aa;A/A) + — / dA — cos(ax A /A) 
cr J u\ a 2 J u\ 

o o 

J A 0J\ 





oo 



^=0 I " n 



2 ^ ^kK/A) ^^ 



o 

?+1 /■ dA uj x - 1 



A 



ge+i(-x\) 



lf(x x ,a/\) 



+ 21+(x x , a/X)li(x x , a/A) - 22? (x A , a/A) 



(53) 



This is the central result of this paper, the explicit expression for the structure factor of the 
OCP, given in terms of an infinite series. To check the internal consistency of the theory, we 
explore how well it satisfies various known sum rules. All these can be conveniently summarized 
in the exact, small k expansion of the structure factor [2], 



'S'exact(^) 



(54) 



where x is the compressibility of the OCP. The first term is the result of the charge neutrality 
and of the Stillinger-Lovett second-moment condition [10], while the second term corresponds to 
the fourth- moment or the "compressibility" sum rule [2] . The inverse compressibility is defined 
thermodynamically in terms of the variation of the pressure, 



(55) 
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with respect to the density, 



= Pttt = 1 + 2 p+jtt + pIt^ = oar, • ( 56 ) 



P+Xp dp + dp + + dp 2 + 36u 



For small couplings this can be expanded to yield 



P = i _ ^?r 3 / 2 + -r 3 - ^r 9 / 2 + — r 6 + o(r 13 / 2 ). (57) 



P+Xp 4 2 2V3 4 

To see if our expression for S(k) is consistent with the sum rules, we expand (53) around a — 0, 
using the asymptotic form of the spherical Bessel function of the first kind, 

. / x (ax) e (axY . nS 

lim jAax) = -\ — J — = K —r-„ 58 

«-o Jn ! (2£ + l)H 1.3.5... (2^ + 1) V ; 



It is evident that, up to order O{k ), only the isotropic (£ = 0) terms of (53) contribute to the 
structure factor. We find 



s-\k) = ( 



k 



P 



P+Xs 



+ 0(k 2 ), 



(59) 



P 



P+Xs 



12 J 



ddJ; 



(2cu A + l)(2cu 2 -cu A + 2) 



Uj\ + UJ X + 1 



= 1 - - 

6 



1 - 2uu + uo' 



TV 



V3 3, (uj 2 + uj + \\ 3^ 



. -In 
2 4 



+ 



tan 



2uj+ 1 

~7T, 



• (60) 



In the low-density limit the inverse compressibility derived from the structure factor can be 
expanded to yield 



A_ = i _ ^! r 3 / 2 + -t 3 



P + Xs I ^ 1 2- 2v ^ +16 1 +U[L h 



(61) 



We see that the structure factor satisfies exactly the charge-neutrality and the second-moment 
conditions, while the compressibility sum rule is satisfied only to order T 9 / 2 . This results from the 
fact that, in order to simplify the calculations, we have neglected the dependence of the cavity 
size and shape on the imposed density variation (A). Clearly, if this was taken into account, 
the theory would be completely internally self consistent. Nevertheless, even at this level of 
approximation, the lack of self consistency is quite small over the full range of relevant coupling 
constants, as can be measured by the inverse compressibilities derived from the thermodynamic 
(Xp) an d the structure factor (xs) routes, Eqs. (56) and (60), respectively; see Fig. (1). The fit 
of the MC data [14] leads to an inverse compressibility (xmc) which is between the two previous 
ones. 



12 





-20 
-40 
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40 80 120 160 

r 

Fig. 1. Comparison between the inverse compressibilities derived from the pressure, Eq. (56), solid 
line (/3/p+xp), and from the structure factor, Eq. (60), long-dashed line ((3/p+xs)- The dashed line 
(P/p + XMc) represents the fit of the MC data over the interval 1 < F < 160 [14]. 

Defining k = kd — aVST, we have explicitly carried out the summation for the first six terms 
of the infinite series in Eq. (53). The results for the structure factor for various values of T, 
obtained without any fitting parameters, are plotted in Figs. (2) to (5). The agreement with the 
MC simulations [15] is quite encouraging. We should note, however, that for higher couplings, 
in the vicinity of the first peak, the series is slowly convergent. This is the reason why we did 
not attempt to carry the calculations for T > 40. 



4 Conclusions 

We have presented the generalized Debye-Hiickel theory of the one-component plasma. The 
linearity of the theory allows for explicit calculations of all the thermodynamic functions, as 
well as the structure factor, which is expressed as an infinite series. The linearity also insures the 
internal consistency of the theory. The agreement with the Monte Carlo simulations, obtained 
without any fitting parameters, is quite good, suggesting that the existence of an effective cavity 
surrounding each ion is theoretically justified. 
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Fig. 2. Structure factor S(k) for T = 2. The solid line is our expression (53) calculated up to £ = 6, 
while the circles represent the MC data [15]. 
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Fig. 3. Structure factor S(k) for T = 6. The solid line is our expression (53) calculated up to £ = 6, 
while the circles represent the MC data [15]. 
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Fig. 5. Structure factor S(k) for T = 40. The solid lines are our expression (53) with the £ = 0, 1, 2, . . . 
up to the £ = 6 terms included in the sum (from top to bottom). The circles represent the MC data 
[15]. We note that for higher values of T, more and more terms will need to be included in order to 
achieve convergence in the vicinity of the first peak. 



A Some useful integrals 



In this appendix we present some integrals which appear along the text. We introduce the 
Kronecker delta, 

4o = ^(27r) 3 5 3 (fc), (A.l) 
15 



where V = (2n) 3 5 3 (0) is the volume of the system and 5 3 (k) is the three-dimensional Dirac 
delta function, 



5 3 (k) = 3 J d 3 r' exp (ife • r') . 



(A.2) 



From the definition (A.2) of the three-dimensional Dirac delta function, it follows directly 



J d 3 r' cos k ■ r' = ^ J dV [exp (ife • r') + exp (-ife • r')\ 

= l - (2tt) 3 [5\k] + 5 3 (-k)} = (2n) 3 5 3 (k) = V5 k0 , 
J d 3 r' cos 2 k-r' = ^J dV (1 + cos 2k ■r') = ^[v+ (27r) 3 5 3 (2fc)] 

= ^(1 + M- 



A simple generalization of (A. 4) leads to 



/dV 



cos 2 • r + 



7T\ 1/ 



! + (-!)%* 



(A.3) 



(A.4) 



(A.5) 



where £ is an integer. 

Expressing 1/r as the inverse Fourier transform 

^/d 3 fclexp(-ifc.r), 



1 

r 



(A.6) 



we have 



/ 



d 3 r dV 



cos k ■ r cos k ■ r' 1 



\r — r 



4tt 2 



f c [3 r j3 r / ^3^ J_ r_-g . ^ r _ r /\j 

J g 2 

x [cos fc • (r + r') + cos k ■ (r — r')\ 

If 1 
g — J d 3 r dV d 3 q — ^ exp [i (anfc - q) ■ r] 



x ex P [i ("2^ + q) ■ r'} 

«2=± 

i(2vr) 4 /d 3 qi ^ 5 3 («ifc-g)^ 3 («2fc + g) 

(27T) 4 

A; 2 
2tt\/ 

"fc 2 " 



5 3 (0) + ^ 3 (2fc) + i<5 3 (-2fc) 
(i + 4o)- 



(A.7) 
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B Variation of the free-energy density 



In this appendix we obtain the variation of the reduced free-energy density, Sf = (35F/V, up 
to quadratic order in the perturbation parameter A. 

The Helmholtz free energy F is written as a functional of the mobile-ion density p+{r). The 
variation 5F is obtained using the functional Taylor series, 



P6F\p+(r)]=PF\p + (r)]-PF [p H 



j3 / P*F 

a r 



5p + (r') 



Sp + (r') 



P+0)=P4 



2 J Sp+(r')Sp+(r") 



8p + {r')8 P+ {r"). 



P+(r)=P+ 



(B.1) 



The linear term can be written as 



/35FU [p+ ( r )} = J dV ^ ^ 5 P+ (r>) = J dV (3p{r') Sp + (r% (B.2) 



p+(r)=p+ 



where p(r) is the chemical potential at the position r, 

5F 



p(r) = 



Sp + (r) 



P+(r)=P+ 



;b. 3 ) 



However, at the thermodynamical equilibrium, the chemical potential of the system is constant 
and is independent of position, 



p(r) = p, 



Vr. 



(B.4) 



Using the imposed variation (7) of the mobile-ion density, 

5p+(r) = p + A cos k ■ r — -p+A [exp (ife ■ r) + exp (— ik ■ r)] , (B.5) 



the linear term can be expressed as 

(35F {1) [p+(r)\ = f3pp + A J dV cosfc • r' = (3pp + VA5 k0 . 



(B.6) 



The quadratic term, 
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5p + (r')5p + (r") 



= ly dVdv 



p+(r)=p+ 



5 (r' - r") 



- C++ (r' - r") 



P+ 



5 P+ (r>)5 P+ (r"), (B.7) 



can be split into two parts, the ideal-gas contribution, 

/^41i [P + (r)\ = ±-J dV dV 5(r' - r") 8p + (r') 8p + (r"), 

and the electrostatic contribution, 

W2L \ P+ {r)} = -\j dV dV C++ (r' - r") 5p + (r') 5p + (r"), 

where C++ (r' — r") is the direct correlation function. 

Using (A. 4) and (B.5), the ideal-gas contribution can be straightforwardly obtained, 



2^/dV[<Wr')f 
^p+A 2 (1 + 4o) • 



ip+A 2 | dV cos 2 fcV 



(B.8) 



(B.9) 



(B.10) 



To evaluate the electrostatic contribution we use (B.5), and express C++ (r' — r") as the inverse 
Fourier transform 



C++(r) = ( J-)' / d 3 fc C++(fc) exp (-ifc ■ r) , 



(B.11) 



leading to 



[p + (r)} = -\p1^ 2 (I)' J dVdV'd 3 g C++(g) 

x ^2 ex P [i ~~ 9) ' r 1 X! ex P [* ( a2 ^ + 9) * r "] 

= -ip 2 + A 2 (2tt) 3 / d 3 q C++(g) ]T 5 3 ( ai fc - g) 5 3 (a 2 fc + g) 

8 J ai,a 2 =± 

= -ip 2 +A 2 (2tt) 3 {C++(fc) [5 3 (0) + 5 3 (2k)] + C++(-fc) [5 3 (0) + 5 3 (-2k)] } 



8 

= -^C++(fc)^A 2 (l + 4o), 



(B.12) 



where we have used the symmetry of the direct correlation function, C++(fc) = C++(— k). 
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Combining all the pieces, the variation of the reduced free-energy density can be written as 



Sf l P+ (r)] = ^SFM [p+{r )} + ^41, [ P+ (r)} + ^5F^ ct \p + (r)] 
= pflp + A5 k0 + J [l - p+C ++ (k)] p+A 2 (1 + 4o) 
= (3pp + A5 k0 + ^S-\k)p + A 2 (1 + 4o) , 



(B.13) 



where S(k) is the structure factor. Note that the linear contribution to the variation has a Kro- 
necker delta (<5fco) factor, which expresses the translational invariance (B.3) of the equilibrium 
chemical potential p of the system. 



C Green's function associated with the induced potential 



In this appendix we obtain the Green's function associated with the differential equation sat- 
isfied by the induced potential 0(r,r'). 

The Green's function G (R, R) associated with (31), where R = r — r', satisfies the homoge- 
neous equation [9] 



Yr - 4© (\R\ -h) G (R, R) = -An5 s {R - R) 



Att 



(C.l) 



with 0(£) the Heaviside step function. 

The general solution of (C.l) can be written as [9] 

oo 

G (R, R) = k d G z ( k d-R) k d #) Pi 



R R 

RRi 



(C.2) 



where Pe(£) denotes a Legendre polynomial. 

( R ■ R'\ 

Replacing (C.2) into (C.l), multiplying both sides by and integrating over the 

y RR J 

angular coordinates 6 and ip, we obtain the equation satisfied by the radial functions Gg (s, s'), 



d 2 2 d 1(1+1) 
+ ~~\ 6 (s — x) — 

as z s as s z 



G e (s,s') = --(2l+l)5(s-s>), (C.3) 



where we have introduced the adimensional variables s = k^R, s' = kbR' and x = n^h. To 
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obtain (C.3) we have used the property of the Dirac delta function, 



5(k b R - k b R') = —5{R - R'), 



(C.4) 



the addition theorem for the Legendre polynomials, 



-R • R'\ 

= Pi [cos 9 cos 9' + sin 9 sin 9' cos(<£> — <//) 

RR J 



(C.5) 



= P e (cos9)P e (cos9') + 2j2 )„ , „{\ PP(cos9) Pp(cos9') cosm(ip - (p'), 



(£ + m)\ 



and their orthogonality, 



d(cos0) P e (cos6) P e (cos6) 



2£ + l 



5ppi. 



(C.6) 



The solutions of (C.3) that are finite for s — > and vanish as s — > oo can be written as 



G £ (*,*') = 



A 12 / + A 13S ^ +1 ), 

^31 ^(s), 
k >W<( S ) + A 33 ^(s), 



for < s < s' < x, 
for < s' < s < x, 
for < s < x < s', 
for < s' < x < s, 
for < x < s' < s, 
for < x < s < s', 



(C.7) 



where the coefficients {A mn } are functions of x and s' to be determined by the boundary 
conditions; ie(s) and ke(s) are the modified spherical Bessel functions of the first and the third 
kinds [16], respectively, 



k(s) = y — i*+i/ 2 (s), 




ke(s) = \l — K i+1/2 (s). 

7TS 



(C.8) 
(C.9) 



Using the symmetry property of the Green's function [17], we can rewrite Eqs. (C.7) as 



G [ P (s, s')=A 1 s e s' e + B 1 



for < s, s' < x, 



(CIO) 
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Gf (s, s')=A 2 s e < k e (s >/ 

.(3) 



for < s< < rr < s>, 



G) j (s,s')=^3^(s)^(s') + -B 3 ^(s<)^(s>), for < a; < s,s', 



(C.11) 
(C.12) 



where s< = min(s, s'), s> = max(s, s'), and the coefficients {A n , B n } depend now only on the 
size of the exclusion hole x. 

The coefficients {B n }, n — 1,3, are obtained by imposing the discontinuity of the derivative of 
Ge(s, s') associated with the Dirac delta function, 



ds 



sGP (a, s>) 



J +e ds 



sG[ n) (s, s') 



21+1 
7/ ' 



(C.13) 



where e is a positive infinitesimal. Using the Wronskian of the modified spherical Bessel functions 
[16], 



W[ke(s),i e (s)] = k e (s)i' e (s) - i e (s)k' e (s) 



(C.14) 



this leads to 



#i = l, 

B 3 = 2£+l. 



(C.15) 
(C.16) 



The coefficients {A n },n = 1,2,3, are obtained by imposing the continuity of Gg(s,s') and of 
its derivative across the spherical surface at s — s' — x, 



G? (s, s>) 



G? (s, s>) 



ds 



Gf(s,s' 

d ^(2) 



s=s'+e=x 



ds 



(2) 

Gf (s, s') 



Gf (s, s') 



(C.17) 
(C.18) 



s=s' +e=x 



and using the following relations of the modified spherical Bessel functions [16] to express ie(x), 
kg,{x) and k'g[x) in terms of ie±i(x) and k£±i(x), 



ar 



■■i e+1 (x)ke(x) + ke +1 (x)ie(x), 



(2£ + 1) ke(x) = xk e +i(x) - xke-i(x), 

{21 + 1) k' t {x) = £k t -i(x) + (£+1) k e+1 (x), 



which yield 



(C.19) 

(C.20) 
(C.21) 



ke-i(x) 



x^k e+1 (xy 



(C.22) 
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21 + 1 

x e+2 k e+l (x) ' 

k+i{x) 



A 3 = (2£+l) 



kt+i(x)' 



(C.23) 
(C.24) 



Therefore, the Green's function G(R,R') is given by the expansion (C.2), with the radial 
functions Gi (s, s') defined by [9] 



G? (s, s>) = 



s e+i x 2i+i ke+v 



x 



G { p{s,s') = {2t + l) S < h{3>] 



x E+2 k e+ i(x) 



for < s, s' < x, 



, for < s< < x < s>, 



Gf (s, s') = (2£ + 1) 



k e+1 (x) 



k e (s)k e (s') + ^(s<)^(s>) 



, for < x < s, s'. 



(C.25) 
(C.26) 
(C.27) 



D The perturbative solution of the induced potential 



In this appendix we obtain the induced potential <p(r,r') recursively, up to order A 2 , at the 
center of the exclusion hole, \r — r'\ = 0. 

Let us obtain the induced potential outside the exclusion hole, \r — r'\ > h, which we will denote 
by (f) > (r,r l ). Clearly, this potential is produced by the charge distribution inside and outside 
the cavity. Let us first calculate the contribution to the potential arising from the charge inside 
the hole, 0> < ' ) (r, r'). Since our final goal is to calculate the potential at the center of the cavity 
to order A 2 , it is sufficient to calculate the induced potential outside the hole to order A, see 
Eq. (33). Using the Green's function G (R, R') derived in appendix C, where R = r — r', we 
find to first order in A, 



1 



<pl < \r,r>) = ^ > < \R + r>,r>) = j 5 J d 3 R g (R r ) G (R, R) 

\R'\<h 



qKB 
D 



J d 3 R {d 3 (R) - p+ [1 + Acosfc • (R + r')]} 



\R'\<h 

E l It - It 
Pe 



RR' 



If 1 3 \ k (K D R) qp+Kv 

k e (K D R) 
x e+2 k i+1 (x) 



A Re 



exp (ifc • r') (21 + 1) 



x 



J d 3 R (k d R'Y exp (ifc • R) P t ( 

\R'\<h ^ 



R ■ R'\ 

RR' ) 



(D.l) 
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recalling that x — K^h. The first term of (D.l) can be simplified using the identities 



k 1 (x) = (l + x)—, 



Ab^d — - 



r 



'1 + x) 



x (1 + x) H — x 3 . 

3 



(D.2) 
(D.3) 



The relation (D.3) is the defining equation for the cavity size, x, Eq. (5). It is important to 
remember that it does not take into account the imposed variation in the ionic density, and as 
result will be responsible for the violation of the compressibility sum rule. 

To simplify the second term of (D.l), we note first that, without loss of generality, we can 
choose the z axis along the k direction, 



cos 9 ■ 
tan<^ : 



k R 

kR ' 
R y 



cos 9' 



k R' 



tan Lp 



kR' ' 

, R' y 



R x R x 

so that the addition theorem for the Legendre polynomials can be written as 



(D.4) 
(D.5) 



( R R'\ 

P e = P e [cos 9 cos 9' + sin 9 sin 9' cos(p - <p')] (D.6) 

y RR J 

= P e (cos9)P e (cos9') + 2j2 j- p PP (cos 9) Pp (cos 9') cos m(ip - ip'). 

m=l \" + m )- 

Performing the integrations over the azimuthal angle tp', only the m = terms survive, 



fi< ( r , r') = -^-xe x k (K B R) - ^ A Re 



(3q 



(3q 



exp (ife • r') ]T (2£ + 1) 



k £ (K D R) 
x e+2 k e+ i(x) 



n i 

x P e (cos9) J dR' R' 2 (k d R') £ j d(cos^)exp (ifci? cos 0') P e (cos9') 



To proceed, we use the plane-wave expansion, 



(D.7) 



exp (ikR'cos9') =Y,( 2 t+ 1 )i i je(kR')Pe(cos9'), 



(D.8) 



where is the spherical Bessel function of the first kind, 



MO = \I^Je+i/2(0- 



(D.9) 
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The integrations over the polar angle 6' and over the radial coordinate Rl can be performed 
using the orthogonality of the Legendre polynomials, 



J d(cos0O^(cos0')^'(cos0') = 2^TT 5 «'> (D.10) 

-i 



and the recursion relation for the spherical Bessel function of the first kind, 
which yields 



[e +2 jt + m]=z e+2 m, (d.ii) 



1 1 OO / 

<g°(r, r') = 4> { <\R + r', r') = —xe x k (K D R) - —A £ (2£ + 1) cos (fc • r' + 1- 

pq otpq V Z 

x J ' +l{( * X h e(K D R)Pt(cosO), (D.12) 

where a = k/ k-q. 

Substituting 0>^(r, r') into the expression for the charge density outside the exclusion hole, 
Eq. (33), we can now calculate the contribution to the potential outside the exclusion hole 
arising from the external charge, 4>> \r,r'). To order A we find 



r,r') = 4> { >\R + r',r') = - J d 3 R' g (R!) G (R, R') 

\R'\>h 

= -^A f d 3 R' + r',r') cos k ■ (R> + r>) 

Att J 



\R'\>h 

oo 

x^Gf^nuR, k d R')P £ 

1=0 

oo 



RR' 



1 OO / \ 

= -— xe x A2 (2£+l) cos (fc • r' + £- ) E £ {k b R, a) P e (cos6), (D.13) 



where the function E^(s, a) is defined by 

,(3) ( 



3<(a,a) = / dsV 2 ^^^^') 

. . oo s 



X X 
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oo 

+ i e (s) J dtek (Oke(OjM). 

s 



(D.14) 



We are now able to find the induced potential inside the exclusion hole, \r — r'\ < h, up to 
order A 2 , 



(r, r') =MR + r', r') = 1 J d 3 R g (R 1 ) G (R, R') 

J d 3 R' {5 3 {R') - p + [1 + A cos k ■ (R + r')} } 



D 



\R'\<h 



R R'\ Kr 



x £ G { P (k d R, k d R') P e ( ± -^ r ) - ^ A / d 3 # 0> (# + r', r') 



£=0 



x cos fe • [R + r') ^ df \k b R, k d R') P e 



\R'\>h 

R R' 

RR' 



where the induced potential outside the exclusion hole, up to order A, is given by 

> (r,r') = 0i <) (r,r') + 0i >) (r,r / ). 



(D.15) 



(D.16) 



Since we need just the mean induced electrostatic potential ip(r') felt by the positive ion fixed at 
r', that is, at the center of the exclusion hole, R = 0, and recalling that Ge(K D R = 0, k, d R') = 
0, W > 0, only the isotropic (£ = 0) terms of (D.15) contribute, 



ib(r') = lim 



4><(R + r',r')- 



DR 



= lim 



J d 3 R'G ( o ) {K B R,K D R')5 3 (R') 



\R'\<h 



DR 



D 



J d 3 R' G ( o ] (0, k d R') [1 + A cos k ■ (R 1 + r')} 



\R'\<h 



^A / d 3 R' G ( q ] (0, KpR^^iR + r^r') cos k-(R' + r') 
An J 



lim 



\Rf\>h 



^jdR'R' 2 G^(0,K B R') 



3 

^A cos fc • r' J dR' R' 2 G { q\0, K D R')j (kR') 



(3q 

4 1 

47r x 2 ki(x) 



A / dlR%(« D iOM# + r V , ) c °sJfe-(- R ' + » ,/ )- (D.17) 



\R'\>h 
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Using the explicit form of Gq~\s, s') and 4>>(r, r'), and performing the angular integrations, we 
obtain 



v ' o 
— A cos k ■ r' < j ds s 



.0 

oo 



1 _ 

s xki(x) 
1 k-i(x) 



+ 



s xki(x 

OO / 

A 2 ^(2£+ l)cos 2 (k-r' + £ 

0—C\ \ 



jo (as) + 



xe 



2x 



1+X 



a(l + x) ^ V ' 2 J ke+i(x) 

,2x oo 



dss k (s)j (as) 



ds s 2 k {s)k e (s)je{as) 



+ ——A 2 £ (2£ + 1) cos 2 (k-r' + £-) ds s 2 k (s)~ e (s, a)j e (as). (D.18) 
1 + x t=o v 2 y ^ 

The first terms of (D.18) can be simplified using (D.2) and (D.3), supplemented by the identities 



k (x) = 



x 



k-i(x) ki(x) — k (x) / x x 
k±(x) ki(x) 1 + 

Furthermore, expressing ie(£) in terms of using the relation [16] 

1 



(D.19) 
(D.20) 
(D.21) 

(D.22) 



and defining the £-th grade polynomial ge(£) associated with the modified spherical Bessel 
function of the third kind ke(£) by the identity [16] 



-^m\\2m)^ ' (D - 23) 



where Y(m) = (m — 1)! is the Euler gamma function, it is possible to express the last integral 
of (D.18), which is two-dimensional, in terms of one-dimensional quadratures, 



J ds s 2 k (s)~ e (s, a) U {as) = e^-l)* 1 ^^^p [^^ a) 



[2 ge+i(x) 

+ If(x, a)li(x, a) - Z?(x, a) | , (D.24) 



where v — ±, 0, are the one-dimensional integrals 
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o 

oo 

J°(x, a) = J ds s~ e g e (s)je(as)2l(s, a) exp [2(x - s)] , 

X 

oo 

J/" (re, a) = j ds s~ e g e (s)j e (as) exp [2(x - s)] . 



(D.25) 
(D.26) 
(D.27) 



We stress that the functions {X^(x,a)} represent one-dimensional quadratures, since the inte- 
grals {Ii(s, a)} can be expressed in explicit form, for all values of £, in terms of trigonometric 
functions and of the sine integral, 



si (t) = f di 



sin£ 



(D.28) 



To illustrate, we give the three first integrals: 



I (s, a) 



Si (as) 



a 



t-/ x 1 cos as ( \ 1 \ . 1 . . 

Jr(5 ' a) = ~a + ^7 + ~ 2^J Sm aS + 2 Sl M ' 

J-( S) a) = -l - 



3 9 3 

+ - — tt^ + — I cos as 



+ 



2a 2 s a 2 s 2 4a 2 s 3 8s, 

3 3 9 3 \ . / 1 3a , . 

2^ + ^-4^ + 8^J SmaS+ fe + T> SlK 



(D.29) 
(D.30) 

(D.31) 



Therefore, the final form for the mean induced electrostatic potential at the center of the 
exclusion hole (in unities of (3q) reads 



(3qip(r') = —-x(x + 2) — A cos k ■ r' 



1 sin ax 



cos ax x , 

+ — — Xo{x,a) 



a 2 a 3 (l + :r) a 2 (l + :r) l + x 

TV 

2 



+ , 1 , A 2 E {2t + 1) cos 2 f fc • r' + ^ +2jm( ^ ) 2+(x, a) 
a(l + x) ^/ ; V 2 7 g e+1 (x) 



x 



1 g e+ i(-x) 

2 &+i(aO 



J/(x,a)l +J/(x,a)J^ (x,a) -l$(x,a) \ . 



(D.32) 
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